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Transport  of  liquid  water  in  a  fuel  cell  gas  diffusion  layer  is  analyzed  using  capillary  network  modeling 
in  which  the  mobility  of  both  liquid  and  gas  phases  is  considered  to  examine  two  distinct  multiphase 
flow  regimes:  displacement  and  co-current  flows.  The  simulations  utilize  a  modified  invasion  percolation 
with  trapping  algorithm,  and  the  capillary  network  consists  of  throats  of  different  radii  to  account  for 
the  local  heterogeneities  of  the  porous  media.  Both  displacement  and  two-mobile  phase  flow  are  solved, 
with  inlet  boundary  condition  for  two-mobile  phase  flow  prescribed  through  a  discrete  sequence  of  alter¬ 
nating  phases  entering  the  network.  For  both  flow  types  (displacement  and  two-mobile  phase),  the  cases 
studied  range  from  capillary  force  controlled,  where  the  maximum  distance  between  two  throats  filled 
consecutively  is  equal  to  the  network  size,  to  viscous  force  controlled,  where  the  maximum  distance 
is  set  so  as  not  to  exceed  some  predefined  value  that  is  less  than  the  network  size.  The  maximum  dis¬ 
tance  determines  the  distribution  of  phases;  phase  entrapment,  percolation,  and  channeling  are  observed 
during  the  spread  of  phases  for  distinct  flow  conditions.  Once  a  distribution  of  phases  is  obtained,  we  cal¬ 
culate  saturation,  relative  permeabilities,  and  the  capillary  pressure  at  the  interface  between  the  phases; 
we  also  determine  the  dependence  of  these  transport  parameters  on  medium  heterogeneity  and  cluster 
size.  Finally,  the  changes  of  relative  permeability  and  capillary  pressure  as  a  function  of  saturation  are 
compared  for  displacement  and  two-mobile  phase  flow. 
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1.  Introduction 

Water  plays  a  critical  role  in  the  operation  of  PEM  fuel  cells.  In 
order  to  achieve  models  that  are  physically  representative  of  the 
multiphase  flow  nature  of  water  transport  in  gas  diffusion  layers 
(GDLs),  two  issues  need  to  be  resolved.  The  first  is  understanding 
of  the  water  flow  mechanisms,  when  liquid  water  displaces  the 
gas  phase,  or  when  gas  and  liquid  phase  flow  co-currently  through 
the  GDL.  The  second  is  the  prediction  of  the  multiphase  transport 
parameters  for  a  given  flow  mechanism.  Two  additional  transport 
processes  need  to  be  kept  in  mind:  water  evaporation  and  capil¬ 
lary  flow,  which  can  greatly  influence  the  liquid  water  flow  pattern 
[1,2].  MRI  experimental  results  [3]  for  sessile  droplets  show  that 
the  evolution  of  imbibed  liquids  can  differ  significantly  depending 
on  the  relative  magnitude  of  evaporation  and  capillary  flow  rates. 

Once  the  water  transport  mechanism  is  defined,  the  liquid  water 
distribution  can  be  obtained  by  solving  the  multiphase  conser¬ 
vation  equations  [4,5].  It  is  generally  accepted  that  the  capillary 
pressure  can  be  expressed  using  the  Leverett  J-function  [6,7]  and 
the  relative  permeability  changes  as  a  power  law  of  the  saturation 
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[8,9],  and  this  is  the  basis  of  most  macroscopic  fuel  cell  models 
accounting  for  liquid  water  transport  in  the  GDL  Although  the 
definitive  form  of  relations  specifically  applicable  to  the  fibrous 
hydrophobic  GDL  media  remains  to  be  established,  some  notable 
progress  in  this  direction  has  been  made  recently  [10,11].  The 
parameters  that  influence  the  relative  permeability  power  law  are 
discussed  and  summarized  by  Valavanides  and  Payatakes  [12], 
where  the  value  of  the  exponent  changes  due  to  the  difference 
in  length  scales  of  the  porous  medium  and  distributed  phase(s) 
[13,14]  that  may  be  caused  by  medium  structure  or  process  dynam¬ 
ics  [15]. 

The  dependencies  between  different  scales  can  be  elucidated  by 
using  capillary  network  models  [16]  in  which  the  porous  medium 
is  represented  as  a  network  of  pores  and  throats.  Three  different 
forces  -  capillary,  viscous  and  gravitational  -  govern  the  multiphase 
flow  [17,18],  and  their  relative  influence  on  the  multiphase  flow  is 
quantified  by  the  capillary  and  Bond  numbers.  For  very  slow  flow 
(small  capillary  number)  and  negligible  gravity,  the  flow  is  dom¬ 
inated  by  capillary  forces  [19].  It  has  been  shown  [20]  that  slow, 
non-wetting  liquid  displacement  flows  follow  the  invasion  perco¬ 
lation  process.  This  percolation  problem  has  been  formulated  by 
Wilkinson  and  Willemsen  [21  ],  whereby  the  interface  advances  at 
each  discrete  step  only  at  one  point  in  which  the  potential  thresh¬ 
old  is  satisfied.  For  a  wetting  liquid,  or  higher  capillary  numbers, 
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dynamic  schemes  [22]  need  to  be  used.  In  this  case,  for  one  discrete 
step,  more  than  one  point  at  the  interface  can  satisfy  the  poten¬ 
tial  threshold,  and  each  pore  can  be  filled/emptied  more  than  once 
during  the  flow  process. 

The  fuel  cell  gas  diffusion  layer  is  an  inherently  heterogeneous 
fibrous  porous  medium  [10]  in  which  liquid  water  is  distributed 
to  form  percolating  flow  paths.  Using  the  percolation  algorithm,  a 
network  structure  which  represents  a  GDL  can  be  defined  to  match 
experimental  capillary  pressure  curves  of  some  actual  GDLs  [23]. 
The  changes  of  the  relative  permeability  and  capillary  pressure  in 
the  limit  of  capillary  flow  as  predicted  from  the  invasion  percola¬ 
tion  are  reported  by  Markicevic  et  al.  [24],  where  this  algorithm 
was  used  to  investigate  directed  water  transport  [25].  Nam  and 
Kaviany  [26]  have  used  network  simulations  to  establish  how  the 
diffusion  coefficient  depends  on  the  saturation,  using  a  cubic  law 
for  the  relative  permeability  and  a  J-function  for  capillary  pressure 
in  the  momentum  transport  equation.  Network  models  reveal  [27] 
that  both  geometrical  and  capillary  properties  influence  the  liq¬ 
uid  water  distribution  in  the  GDL.  Experimental  evidence  [28-30] 
shows  that  liquid  water  clusters  throughout  the  GDL;  thus,  water 
emerges  in  specific  points  in  the  channel  area  rather  than  flow  con¬ 
tinuously  across  the  overall  available  interface  between  the  GDL 
and  the  channel.  Based  on  the  cluster  size,  the  correlation  length  of 
the  liquid  water  flow  can  be  calculated,  where  recently  it  has  been 
shown  [31]  that  the  GDL  thickness  is  smaller  than  the  correlation 
length  which  limits  the  use  of  continuum  models. 

In  this  study,  co-current  flow  of  two  mobile  phases  is  investi¬ 
gated  and  the  impact  of  variations  of  flow  rates  of  both  phases  at 
the  inlet  of  the  porous  medium,  of  the  maximum  allowable  clus¬ 
ter  size  (correlation  length),  and  of  the  heterogeneity  of  the  porous 
medium  are  examined.  For  the  limiting  case  when  one  of  the  phases 
at  the  inlet  boundary  is  dominant,  the  co-current  flow  reverts  to  dis¬ 
placement  flow.  For  the  spread  of  both  phases,  invasion  percolation 
with  trapping  is  used,  whereby  both  phases  form  static  structures 
and  momentum  transfer  between  the  phases  is  negligible.  The  trap¬ 
ping  rule  is  altered  so  that  it  can  account  for  the  presence  of  two 
mobile  phases.  Furthermore,  the  invasion  percolation  algorithm  is 
modified  for  variable  correlation  length,  in  which  the  width  of  the 
advancing  interface  cannot  exceed  a  predefined  correlation  length. 
The  invasion  is  stopped  at  the  termination  point,  when  all  throats 
at  the  outlet  boundary  are  invaded  by  one  of  the  phases  transported 
from  the  network  inlet.  Based  on  the  resulting  patterns,  the  rela¬ 
tive  permeability  of  each  phase  as  well  as  the  capillary  pressure  at 
the  interface  are  calculated  and  the  impact  of  saturation  and  other 
investigated  parameters  is  analyzed. 


2.  Model  system 

A  discrete  capillary  network  model  for  two-phase  flow  in  which 
both  phases  are  mobile  is  developed.  The  network  consists  of 
throats  and  pores,  where  the  network  coordination  number  rep¬ 
resents  the  number  of  pores  connected  to  one  pore.  For  regular, 
square  networks,  the  network  coordination  number  is  constant  and 
equal  to  four.  A  schematic  of  the  network  and  definitions  of  the 
parameters  are  shown  in  Fig.  1.  In  the  network,  each  throat  can 
only  be  occupied  by  one  phase.  Both  phases  enter  the  network  at 
the  inlet  boundary  and  in  order  to  vary  their  flow  rates,  the  number 
of  invading  steps  of  the  first  (ni )  and  second  phase  (n2 )  is  prescribed 
(e.g.  when  rq  =1  and  n2  =  2,  for  each  throat  invaded  by  the  first 
phase,  two  throats  are  invaded  by  second  phase).  Hence,  the  occur¬ 
rence  frequency  is  defined  as  a  ratio  of  invading  steps  (v  =  ni/n2), 
and  either  phase  can  be  dominant  compared  to  the  other.  The  pro¬ 
cess  repeats  alternatively  until  there  are  no  more  available  throats 
to  occupy;  this  corresponds  to  the  termination  point.  Once  the 
flow  paths  of  the  phases  are  initiated,  they  branch  into  the  porous 
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Fig.  1.  Network  model  of  a  porous  medium  consisting  of  pores  and  throats.  The 
flow  parameters  are  shown:  occurrence  frequency  of  the  phases  at  the  network 
inlet  (v  =  ni/ri2),  medium  heterogeneity  ( X  =  rmaxlrmin  -  1),  and  the  length  scale  rep¬ 
resenting  the  cluster  size  (£  =  /s/L). 


medium  and  form  distinct  flow  patterns.  The  irregularity  of  the 
pattern  of  each  phase  is  caused  by  local  heterogeneities  that  pro¬ 
duce  local  variations  of  permeability  (I<)  and  capillary  pressure  (pc), 
which  in  turn  result  in  local  flow  rate  variations  for  each  phase.  The 
heterogeneity  is  determined  by  the  prescription  of  the  throat  radii 
(rt)  using  a  predefined  distribution  function.  For  a  uniform  distri¬ 
bution  of  throat  radii  (rt)  with  a  minimum  (rmln)  and  maximum 
{rmax),  Tm/n  <rt<rmax,  the  heterogeneity  parameter  can  be  defined 
as  x  =  rmaxl^mm  ~  1 .  and  in  the  limit  of  a  homogeneous  sample  we 
have  x  =  0.  It  should  be  noted  that  a  pore  network  does  not  provide 
a  direct  geometric  representation  of  a  physical  porous  medium. 
Physical  consistency  is  ensured  by  selecting  a  throat  radii  distri¬ 
bution  function  such  that  I<  and  pc  in  the  numerical  network  are 
physically  representative  of  values  measured  in  the  actual  porous 
medium  of  interest. 

Viscous  forces  influence  the  length  scale  for  which  the  parts 
of  one  phase  are  fully  surrounded  by  the  other.  The  size  of  the 
entrapped  phase  (cluster)  is  a  measure  of  the  viscous  correlation 
length.  Still,  the  capillary  force  controls  the  phase  displacement  on 
the  local  level  [32].  Therefore,  a  length  scale  (Zs)  is  defined  that  is 
smaller  or  equal  to  the  domain  length  (L)  in  the  network  where 
the  invasion  of  phases  occurs  (note  that  ls  differs  from  the  corre¬ 
lation  length).  Once  some  fraction  of  the  throats  is  invaded  over  a 
length  ls,  the  throats  in  the  next  layer  of  the  network  of  length  ls 
are  allowed  to  be  filled.  This  is  repeated  until  all  layers  of  length  ls 
in  the  network  are  filled.  Hence,  for  a  length  scale  ls  which  is  equal 
to  the  length  of  the  domain  (L),  the  process  is  governed  by  capillary 
forces  and  corresponds  to  a  small  capillary  number.  As  the  length 
scale  ls  decreases,  the  filling  is  increasingly  controlled  by  viscous 
forces,  which  coincides  with  a  larger  capillary  number.  To  vary  the 
length  scale  ls,  the  process  extent  parameter  (£  =  Zs/L)  is  defined  and 
used  (see  Fig.  1). 

Satisfying  the  condition  that  the  spread  of  phases  is  governed 
by  capillary  forces  at  the  local  level,  the  invasion  percolation  with 
trapping  algorithm  [21]  is  used  for  both  phases.  Each  phase  occu¬ 
pies  different  throats  in  the  network  and  there  is  no  momentum 
transport  between  them.  At  the  onset,  the  spread  of  both  phases 
starts  into  a  network  completely  filled  by  the  first  phase  (originally 
present).  This  first  phase,  together  with  the  other  (second)  phase, 
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enters  the  network.  As  two  mobile  phases  are  present,  different 
types  of  clusters  can  be  formed,  and  the  trapping  rules  must  be 
modified.  The  set  of  rules  for  how  the  phases  are  allowed  to  spread 
within  the  network  are: 


(i)  for  each  of  the  two  phases,  the  carrying  backbone  is  defined  as 
a  flow  path  that  spans  the  network  inlet  and  outlet, 

(ii)  the  paths  of  the  different  phases  can  intersect  within  the  net¬ 
work, 

(iii)  only  clusters  of  the  originally  present  phase  in  the  network 
can  be  formed,  and  they  are  surrounded  by  either  of  the  two 
phases’  preferential  flow  paths. 


The  clusters  which  are  surrounded  by  the  originally  present  (first) 
phase  flow  paths  contribute  to  momentum  transport  (they  are  part 
of  the  carrying  backbone  of  the  originally  present  phase)  and  they 
are  mobile.  In  this  case,  mobile  means  that  the  fluid  phase  is  mobile, 
but  the  corresponding  clusters  stay  in  the  same  position.  On  the 
other  hand,  clusters  of  the  first  phase  which  are  surrounded  by  flow 
paths  of  the  other  (second)  phase  do  not  contribute  to  momentum 
transport,  and  in  these  clusters,  the  first  phase  is  immobile.  Clusters 
of  the  second  phase  cannot  form;  however,  there  are  still  parts  of 
the  second  phase  that  do  not  belong  to  the  second  phase  backbone 
and  do  not  contribute  to  momentum  transport.  Finally,  throats  that 
belong  to  the  clusters  which  are  formed  by  first  phase  flow  paths 
(mobile  clusters)  can  be  invaded  by  the  second  phase. 

Within  the  network,  there  are  parts  of  each  phase  that  belong  to 
each  phase  carrying  backbone  (mobile)  and  parts  that  are  immo¬ 
bile.  Thus,  the  content  of  each  phase  is  given  as  a  phase  saturation 
Si  and  s2  =  1  -  Si,  and  the  immobile  parts  are  defined  as  immobile 
saturations  slm  i  and  sim >2.  Since  each  throat  is  occupied  by  only  one 
phase,  the  saturation  is  a  ratio  of  the  volume  of  throats  filled  with 
a  particular  phase  (Vtl)  to  the  volume  of  all  throats  {Vtau): 


vt,i 

Vt,all  ’ 


(1) 


where  the  index  i  refers  to  the  phases  (i  =  1 ,  2),  and  the  index  j 
is  a  summation  index.  Similarly,  the  immobile  saturation  can  be 
obtained  by  adding  the  volumes  of  throats  filled  by  immobile  phase. 

In  single-phase  flow,  the  network  as  depicted  in  Fig.  1  shows 
the  flow  resistance,  which  is  defined  as  a  single  phase  permeabil¬ 
ity  ( Ksp )  or  network  permeability,  and  is  calculated  from  Darcy’s 
law.  This  permeability  is  a  function  of  porosity,  I<sp  =  I<sp (0).  With 
two  phases  present,  the  flow  resistance  of  each  phase  becomes 
larger  compared  to  single-phase  flow.  For  multiphase  flow,  it  is 
assumed  again  that  the  phase  velocity  (Uj)  and  the  phase  pressure 
gradient  ( vp* )  are  related  linearly  according  to  Darcy’s  law.  The  pro¬ 
portionality  constant  is  phase  specific  and  depends  on  the  phase 
permeability  (Kj)  and  the  phase  viscosity  (/xz)  (i  =  1  or  2).  Thus,  for 
each  phase  [33]: 


At  the  interface  between  two  phases,  a  potential  difference 
across  the  interface  can  be  defined.  In  order  for  one  phase  to  spread 
into  another,  it  is  not  sufficient  to  have  a  potential  difference  at  the 
interface,  but  the  difference  must  also  be  larger  than  the  threshold 
potential.  Furthermore,  due  to  heterogeneities  of  the  local  medium, 
the  threshold  potential  varies  both  along  the  interface,  as  well  as 
with  the  motion  of  the  interface  during  the  invasion.  In  the  net¬ 
work,  the  potential  difference  is  given  by  the  pressure  difference 
between  two  connected  pores.  The  threshold  potential  difference 
is  the  capillary  pressure  of  the  connecting  throat  (pc),  calculated 
from  the  Young-Laplace  equation: 


where  o  is  the  surface  tension  of  the  originally  present  phase,  and  rt 
is  the  radius  of  the  throat.  This  displacement  criterion  requires  the 
solution  of  the  pressure  field  during  invasion.  Flowever,  in  the  limit 
of  very  slow  processes  (capillary  force  dominated)  this  criterion 
turns  into  a  simpler  rule:  in  a  drainage-like  process,  the  throat  with 
the  largest  radius  will  be  invaded.  For  faster  flows  (viscous  force 
dominated)  the  same  rule  (largest  throat  radius)  for  the  spread  of 
phases  can  be  exploited,  but  the  maximum  cluster  size  decreases 
as  the  phase  velocity  increases  (the  cluster  size  is  limited  using  ls). 

3.  Numerical  procedure 

A  two-dimensional,  regular  square  network  of  size  (nL  x  nL)  is 
defined  using  throats  with  circular  cross-sections  and  four  throats 
connected  to  each  pore.  Two  phases,  one  of  which  is  the  same  as 
the  initial  phase,  enter  the  network  at  the  inlet  boundary.  The  out¬ 
let  is  on  the  opposite  side  of  the  network.  There  is  no  flow  across 
the  remaining  two  network  boundaries.  Phases  enter  the  network 
alternately,  and  the  number  of  consecutive  steps  of  each  phase  is 
n\  and  n2,  respectively.  If  one  phase  is  dominant,  it  might  block 
the  other  phase,  preventing  it  from  reaching  the  network  outlet.  In 
this  case,  a  deficient  phase  becomes  a  large  stationary  cluster  in  the 
network. 

Once  there  are  no  available  throats  to  be  filled,  the  pressure  for 
each  phase  within  the  network  is  calculated.  It  is  assumed  here  that 
the  phases  are  in  equilibrium.  Thus,  there  is  no  transport  across 
the  interface  (the  interface  is  stationary),  and  each  phase  flows 
only  from  the  network  inlet  to  the  outlet.  Only  the  phase  carry¬ 
ing  backbone  contributes  to  the  momentum  transport  and  to  the 
phase  relative  permeability.  If  one  phase  does  not  reach  the  outlet, 
this  phase  is  regarded  as  a  cluster,  and  a  solution  of  the  pressure 
is  not  required.  For  the  phase(s)  flowing  through  the  network,  the 
pressure  solution  is  found  by  applying  a  mass  balance  at  each  pore 
within  the  network.  The  balance  for  one  pore  (z)  is  given  in  the 
form  of  the  flow  rates  through  all  connected  throats  (q^)  and  has 
the  form: 


u‘  =  “7 7Vp”  (2) 

where  the  phase  permeability  (JQ)  is  a  function  of  porosity  (0)  and 
saturation  (s),  K\  =  Ki(</>,  s).  In  order  to  separate  the  influence  of  the 
saturation,  the  relative  permeability  (/<ri)  of  each  phase  (z  =  1,  2) 
is  defined  as  a  ratio  of  the  phase  permeability  to  the  single  phase 
permeability: 


^r,i 


I<sP  Q ' 


(3) 


As  shown  in  the  above  relation,  given  an  applied  pressure  drop 
in  the  network  (A p/L),  the  relative  permeability  can  be  obtained 
as  a  ratio  of  the  phase  flow  rate  (Qj)  to  the  single-phase  flow  rate 
through  the  network  (Q). 


5>J  =  °.  c  =  4,  (5) 

j=  1 

where  c  is  the  flow  coordination  number.  In  single-phase  flow,  the 
flow  and  network  coordination  numbers  are  equal.  Flowever,  in 
multiphase  flow,  the  flow  coordination  number  is  defined  for  pores 
that  belong  to  the  phase  carrying  backbone.  In  general,  c  varies 
throughout  the  backbone,  where  for  any  pore,  one  or  more  throats 
can  be  occupied  with  another  phase.  The  flow  rate  (q^)  through 
each  throat  (j  =  1, . . .,  c)  is  calculated  from  its  conductance  ( gj )  and 
the  pressure  difference  (A  p)  at  the  starting  and  ending  pore  of  the 
throat  as: 

q,,i=gjAP-  (6) 
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permeability  ( I<sp )  can  be  computed  from  Darcy’s  law.  Similarly, 
both  phase  permeabilities  (Kj)  are  obtained  from  known  phase  flow 
rates  (Qj).  The  relative  permeabilities  are  then  calculated  from  Eq. 
(3).  The  flow  through  the  network  is  two-dimensional,  whereas 
all  permeabilities  are  calculated  for,  in-principle,  one-dimensional 
flow.  Therefore,  for  a  fluid  of  viscosity  (/x): 


k  =  h^2~, 
^  A  Ap 


(11) 


where  I<  and  Q.  are  used  for  either  single  phase  (sp)  or  multiphase 
flow  (i).  The  pressure  drop  in  the  network  is  (pin  -  p0Ut)/L  and  u  =  QJA 
is  the  superficial  velocity.  Finally,  the  capillary  pressure  (pc)  is  found 
by  averaging  the  local  capillary  pressures  (Eq.  (3))  over  the  entire 
static  interface  in  the  spread  termination  point. 


(ij-V 

Fig.  2.  Network  node  with  one  pore  and  connecting  throats. 

For  a  circular  throat  of  radius  (rtJ  )  and  length  (Zt),  the  exact  value  of 
the  conductance  is  obtained  from  the  Poiseuille  solution: 


where  /x  is  the  viscosity  of  the  phase  for  which  the  pressure  solution 
is  solved  (/xi  ^  /x2). 

Applying  the  balance  equation  to  each  pore,  Eq.  (5)  can  be 
assembled  over  the  whole  network  (single  phase  flow)  or  phase 
carrying  backbone.  Fig.  2  shows  a  pore  (i,  j)  and  four  connecting 
throats:  (z-l,j),  (z  +  l,j),  (i, j  —  1 )  and  (z,j  + 1).  Both  horizontal  (H) 
and  vertical  (V)  throats  can  be  previous  (pre)  or  next  (nex)  to  the  pore 
(z,  j).  Since  the  direction  of  the  flow  for  each  throat  is  not  known 
a  priori ,  all  (qj)  are  defined  as  the  flow  into  pore  (z,  j).  Using  this 
convention  and  Eq.  (5),  the  mass  balance  for  pore  (z,  j),  obeys: 

SpreP- 10  +  SnexP- 1-10  +  QpreP- 01  T  SnexP- 1-01 

— (Spre  3" Snex  Qpre  T  Snex^Poo  —  0?  (8) 

where  the  subscripts  for  pressure  (p)  and  conductance  (g)  are  in 
reference  to  Fig.  2.  In  multiphase  flow,  the  throats  that  do  not  belong 
to  the  carrying  backbone  of  the  current  phase  are  omitted,  as  they 
do  not  contribute  to  the  balance  of  the  current  phase.  Applying  Eq. 
(8)  to  each  pore  yields  a  linear  system  of  algebraic  equations  for  the 
pressure  (p)  over  the  network/carrying  backbone.  Hence,  one  can 
write: 

A;P,  =  b„  (9) 

with  index  z  referring  to  the  single  phase  (i  =  sp)  or  multiphase  flow 
(phases,  z  =  1, 2).  The  matrix  (A/)  is  obtained  from  the  conductances 
of  throats,  whereas  the  force  vector  (b;)  depends  on  both  throat 
conductances  and  the  external  pressure  at  the  inlet  and  outlet  of 
the  network. 

Eq.  (9)  is  solved  once  for  single-phase  flow  over  the  whole  net¬ 
work  and  twice  for  the  two  phases  over  corresponding  carrying 
backbones.  Having  determined  the  pressure  distribution,  the  flow 
rate  is  calculated  by  summing  the  flow  rates  through  the  throats  at 
the  inlet  or  the  outlet  boundary.  Once  an  equilibrium  distribution 
of  phases  in  the  network  is  achieved,  mass  conservation  requires 
that  the  flow  rate  of  each  phase  into  the  network  equals  the  flow 
rate  out  of  the  network.  For  phase  z,  the  flow  rate  at  the  inlet  is: 

Qjnlet  =  ^  ^ Sm(Pinl  ~  Pm),  (10) 

m 

where  m  is  a  summation  index  for  all  throats  occupied  by  phase  z. 
From  the  single-phase  flow  rate  (QsP),  the  network  or  single-phase 


4.  Results  and  discussion 

The  numerical  results  are  obtained  for  a  regular  square  network 
(nLxnL)  of  size  nL  =  60.  The  cylindrical  throats  have  dimensions 
(rt,  lt),  where  the  throat  length  /t  =  2xl0-4m.  The  heterogene¬ 
ity  parameter  (x),  the  occurrence  frequency  (v  =  rzi/rz2),  and  the 
process  extent  (f  =  /s/L)  are  altered  in  order  to  vary  the  regimes 
and  obtain  a  broad  range  of  phase  distributions  within  the  net¬ 
work.  Local  heterogeneities  are  included  by  setting  the  throat 
radius  as  a  random  variable.  In  this  study,  the  radii  are  uni¬ 
formly  distributed  in  the  range  (rmint  rmax )  with  an  average  radius 
rav  =  4  x  10-4  m.  The  minimum  (rmin)  and  maximum  (rmax)  radii  are 
set  to  (rmintrmax)  x  10-4m  =  {(3.5,  4.5),  (2.0,  6.0),  (0.5,  7.5)},  with  a 
corresponding  heterogeneity  parameter  (x  =  rmaxlrmin  -  1 )  equal  to 
X  =  {0.3,  2,  14}  ranging  from  almost  homogeneous  to  highly  het¬ 
erogeneous  media.  Either  phase  -  water  or  air  -  can  be  in  excess  or 
deficient  compared  to  the  other  phase,  and  therefore,  seven  distinct 
values  of  v  are  used,  v  =  (0.1 , 0.25, 0.5, 1 , 2, 4, 1 0}.  The  maximum  dis¬ 
tance  (ls)  between  two  throats  that  are  filled  successively  is  initially 
set  equal  to  the  network  length  (L)  which  corresponds  to  capil¬ 
lary  force  controlled  regime.  This  distance  is  gradually  decreased 
to  Zs/L  =  1/20  for  fast  water  transport  corresponding  to  the  viscous 
force  controlled  regime;  thus  f  =  {1, 1/3, 1/10, 1/20}.  Overall  (v,  x> 
£)  =  (7  x  3  x  4)  =  84  different  combinations  of  parameters  are  inves¬ 
tigated.  Finally,  for  each  combination  (v,  x>  ?)>  the  throat  radii  (rt) 
are  randomly  generated  N  =  20  times,  and  the  results  of  all  realiza¬ 
tions  are  used  to  obtain  statistical  averages. 

Once  the  distribution  of  the  phases  is  obtained,  the  saturation  is 
calculated  from  Eq.  ( 1 )  by  summing  the  volumes  of  throats  occupied 
by  each  phase.  Similarly,  the  immobile  saturation  of  each  phase  is 
found,  but  only  for  the  throats  in  which  flow  rate  is  equal  to  zero. 
The  capillary  pressure  is  calculated  at  the  phase  interfaces.  As  the 
spread  of  both  phases  is  locally  governed  by  capillary  forces,  pc  is 
found  by  averaging  the  capillary  condition  (Eq.  (4))  at  the  interfaces. 
The  single-phase,  and  the  phase  permeabilities  are  found  from  the 
phase  flow  rates  across  the  boundary  (Eq.  ( 1 0))  and  Darcy’s  law  (Eq. 
(11)). 

4.1.  Initial  analysis 

A  typical  solution  containing  the  distribution  of  the  two  phases 
is  shown  in  Fig.  3(a);  the  originally  present  phase  (first  phase)  is 
shown  with  gray  lines,  whereas  the  second  phase  is  depicted  by 
black  lines.  The  carrying  backbones  of  both  phases  are  shown  in 
Fig.  3(b).  In  this  case,  the  numbers  of  invading  steps  are  equal  (n\  =  1 
and  rz2  =  1 )  and  occurrence  frequency,  v  is  equal  to  one.  The  length 
scale  (ls )  is  equal  to  the  length  of  domain  (L),  and  f  =  1S/L  =  1 .  As  inva¬ 
sion  percolation  with  trapping  is  used  for  the  distribution  of  phases, 
the  heterogeneity  parameter  (x)  does  not  influence  the  spread  of 
phases.  The  next  throat  to  be  invaded  has  to  be  the  largest  among 
all  available  throats.  As  can  be  observed  from  Fig.  3,  the  distribu- 


B.  Markicevic,  N.  Djilali  /  Journal  of  Power  Sources  196(2011)  2725-2734 


2729 


Fig.  3.  Two  mobile  phases  invasion  of  a  capillary  network:  (a)  distribution  of  phases 
in  the  termination  point;  (b)  phase  carrying  backbones.  The  originally  present  (first) 
phase  and  the  second  phase  are  represented  by  gray  and  black  lines,  respectively. 
The  inlet  is  at  the  top  of  network. 


a 


Fig.  4.  The  invasion  percolation  with  defined  cluster  size  and  trapping  for  the  maxi¬ 
mum  cluster  size  equal  to  one  third  of  the  network  size:  (a)  after  filling  of  first  third; 
(b)  after  filling  of  second  third. 


tion  of  the  two  phases  is  asymmetric.  Although  v  =  1  and  one  would 
expect  equal  saturations  of  each  phase  (si  =s2  =  0.5),  the  saturation 
of  the  originally  present  phase  is  larger  (in  this  realization  si  «  0.55 
and  s2^0.45).  This  difference  is  caused  by  the  formation  of  clus¬ 
ters,  as  the  second  phase  cannot  spread  into  the  immobile  clusters 
of  the  first  phase  and  one  obtains  Si  >v/(v  +  l)  and  s2<l/(v  +  l). 
The  throats  associated  with  the  immobile  first  phase  consist  of: 

(a)  throats  surrounded  by  the  flow  path  of  the  second  phase,  and 

(b)  paths  of  the  first  phase  that  are  not  part  of  the  carrying  back¬ 
bone.  On  the  other  hand,  the  throats  associated  with  the  immobile 
second  phase  belong  to  the  second  phase  paths  only,  but  not  to 
the  second  phase  carrying  backbones;  hence  sim  l  >sim>2.  This  is 
related  to  the  momentum  transport  associated  with  the  first  phase 
(mobile)  clusters.  These  are  formed  as  parts  of  the  first  phase,  which 
are  surrounded  by  preferential  flow  paths  of  the  same  phase  (see 
Fig.  3(b)).  The  mobile  regions  have  low  flow  resistances  which  cause 
two  relative  permeabilities  to  differ  (kr  i  >  kr<2 ). 

Once  the  capillary  number  increases,  the  fluid  spread  is  more 
localized  and  only  throats  which  are  within  a  distance  of  ls  may  be 
filled  in  two  consecutive  steps.  This  stabilizes  the  fluid  front  and 
the  flow  becomes  dominated  by  viscous  forces.  As  a  result  of  flow 
stabilization,  smaller  clusters  are  formed.  Fig.  4  shows  the  filling 
of  the  network  for  f  =  /s/L  =  l/3  after  the  first  (Fig.  4(a)),  and  the 
second  (Fig.  4(b))  layer  of  the  network  is  filled.  However,  changes 


in  ls  do  not  only  alter  the  cluster  size,  but  also  cause  segregation 
of  the  phases.  This  segregation  allows  two  processes  to  occur:  (i) 
channeling  of  phases  and  (ii)  no-percolation  of  the  deficient  phase. 
The  first  process  is  promoted  by  decreasing  ls  and  the  latter  one  is 
observed  for  large  values  of  ls.  In  order  to  quantify  the  percolation 
of  the  phase,  the  dimensionless  fluid  front  position  is  defined  as  a 
ratio  of  the  furthest  fluid  point  from  the  inlet  ( lfront )  to  the  length 
of  the  network  (L): 

(12) 

Thus,  the  phase  percolation  occurs  once  =  1,  and  the  condition  of 
no-percolation  corresponds  to  A;  <  1.  It  can  readily  be  seen  that  for 
the  phase  that  does  not  percolate,  the  phase  saturation  (Sj)  is  equal 
to  the  immobile  phase  saturation  (sim>1),  and  the  relative  phase  per¬ 
meability  is  equal  to  zero. 

In  an  ideal  case,  for  multiphase  parallel  flow  in  which  clusters 
do  not  form,  saturation  is  proportional  to  the  occurrence  frequency 
(v  =  rii/n2),  and  can  be  expressed  as  Si  =  v/(v  +  l)  and  s2  =  l/(v  +  l). 
This  condition  is  satisfied  for  very  low  (f ),  which  corresponds  to 
the  flow  dominated  by  viscous  forces.  Otherwise,  as  £  increases, 
the  influence  of  capillarity  becomes  significant,  and  we  obtain 
Si  >  v/(v  + 1)  and  s2  <  l/(v  + 1).  This  is  shown  in  Fig.  5(a),  where  the 
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Fig.  5.  Estimation  of  the  second  phase  maximum  saturation  (smax, 2)  for  different 
heterogeneity  and  process  extent  (x,  £):  (a)  extrapolation  for  x  =  0.3  and  all  £  inves¬ 
tigated  with  dash  line  limit  s2  =  1/(1  +  v),  and  (b)  power  law  dependence  of  smax>2  as 
a  function  of  the  cluster  size. 


plots  s2  vs.  1/(1  +  v)  for  constant  heterogeneity  x  =  2,  and  for  var¬ 
ious  £  =  /s/L  are  depicted.  The  dashed  line  represents  S2  =  l/(1  +v), 
and  it  can  be  observed  that  the  actual  value  of  s2  is  lower  than  the 
limiting  value  1/(1  +v).  However,  as  the  heterogeneity  increases, 
s2  approaches  1/(1  +  v)  as  the  throats  that  are  filled  are  larger  and 
contribute  the  most  to  the  medium  volume.  The  second  phase  sat¬ 
uration  is  extrapolated  as  v  tends  to  zero  and  1/(1  +  v)  goes  to  the 
intrinsic  saturation  value.  This  saturation  is  referred  to  as  the  max¬ 
imum  achievable  saturation  (smax>2)  of  the  second  phase;  smax>2 
depends  on  £  and  x,  and  the  results  are  depicted  in  Fig.  5(b).  In 
Fig.  5(b),  the  symbols  represent  the  numerical  results,  and  the  solid 
lines  represent  the  power  law  fits: 

Smax, 2  —  01 .  (13) 

The  results  yield  exponents  in  the  range  a  ^0.1 7-0.2,  with 
a  decreasing  with  x*  The  existence  of  the  maximum  sat¬ 
uration  of  the  second  phase  (smax2)  suggests  that  there  is 
an  additional  property  of  two-phase  flow:  the  content  of  the  origi¬ 
nally  present  phase  that  cannot  be  displaced  from  the  network  even 
when  the  second  phase  is  dominant  (in  this  case  v  is  very  small). 
This  retention  of  the  first  phase  is  caused  by  its  immobile  clusters, 
around  which  the  second  phase  flows.  The  amount  of  the  phase  that 
cannot  be  displaced  can  be  quantified  as  an  irreducible  saturation 
(Sir ,1 )  of  the  first  phase.  The  value  of  sir  i  is  related  to  smax>2,  and  the 
irreducible  saturation  of  the  originally  present  phase  (first  phase) 
is: 

$ir,  1  —  1  —  Smax, 2?  (14) 

where  smax>2  is  estimated  from  the  dependence  s2  and  v  as  shown 
in  Fig.  5(a). 


4.2.  Single  (kr,i~s-Pc)  and  (Xt —s)  curves 

In  order  to  assess  how  the  relative  permeabilities  (krl),  capil¬ 
lary  pressure  (pc),  and  fluid  fronts  (A2)  of  each  phase  change  as  a 
function  of  saturation  (s2),  simulations  were  performed  with  con¬ 
stant  values  for  x  and  £,  and  by  varying  v  =  n\\n2.  For  each  value 
of  v,  the  distribution  of  throat  radii  (rt)  is  generated  randomly 
(N= 20)  to  obtain  stochastic  averages.  Hence,  140  discrete  points  (v, 
N)  =  ( 7  x  20)  were  used  to  generate  (kri- s2—  pc)  and  (A.,—. s2)  curves. 
The  relative  permeability  kr  1  is  equal  to  the  permeability  ratio 
K\  lKsp,  where  the  phase  permeability  (JCi)  is  found  using  the  gen¬ 
eralized  Darcy’s  law  (Eq.  (11))  from  the  known  pressure  drop,  and 
calculating  the  phase  flow  rate  at  the  network  flow  boundary  (inlet 
or  outlet,  Eq.  (10)).  The  second  phase  relative  permeability  (kr>2) 
is  determined  in  the  same  manner.  The  saturation  of  the  origi¬ 
nally  present  phase  (first  phase,  Si ),  and  the  immobile  saturations 
(sim  l )  and  (Sjm>2)  are  computed  from  the  whole  network.  The  cap¬ 
illary  pressure  is  averaged  at  the  interface  between  the  second 
phase  carrying  backbone  and  the  remaining  part  of  the  network. 
The  results  for  kr>\ ,  pc,  and  for  £  =  1  and  x  =  14  are  shown  in 
Fig.  6(a)-(c).  The  circles  represent  the  discrete  realizations,  and  the 
bold  squares  show  the  averages  for  each  value  of  v.  The  standard 
deviations  are  also  shown  (error  bars)  in  the  figure,  with  the  hor¬ 
izontal  error  bars  representing  the  saturation  deviation.  It  can  be 
seen  from  Fig.  6(a)  that  for  some  random  network  corresponding 
to  a  small  value  of  v  (the  second  phase  is  dominant),  the  relative 
permeability  kr  \  is  equal  to  zero  ( i.e .  there  is  no  percolation  of  the 
first  phase).  As  v  increases,  the  first  phase  percolates,  and  Si  and 
kr ti  increase.  Therefore,  as  long  as  kr 1  =  0,  the  fluid  front  (k\ )  is  less 
than  one.  The  effect  of  Si  on  the  fluid  front  is  shown  in  Fig.  6(c).  As 
Si  increases,  the  first  phase  protrudes  more  into  the  network  and 
\\  increases  and  eventually  becomes  equal  to  one.  Remarkably, 
Fig.  6(b)  exhibits  a  relatively  constant  capillary  pressure  irrespec¬ 
tive  of  saturation.  This  implies  that  for  specified  values  of  x  and  £, 
the  distribution  of  throat  sizes  at  the  interface  remains  constant 
regardless  of  changes  of  the  flow  patterns  and  interface  shapes 
with  si . 

A  power  law  is  usually  used  to  correlate  the  relative  perme¬ 
ability  (krj)  and  the  saturation  (s2)  in  a  form  kri  =  Bs^.  Such  a 
comparison  for  a  heterogeneity  (x  =  14)  and  all  four  investigated 
£  is  given  in  Fig.  7(a)  and  (b).  In  Fig.  7,  the  symbols  represent  aver¬ 
aged  numerical  results,  and  the  solid  £  =  1/20  and  dash-dot  £  =  1 
lines  are  the  lower  and  upper  power  law  approximations.  For  the 
originally  present  phase  (first  phase),  the  value  of  the  exponent 
(P)  decreases  from  ft  =  2.8  to  ft  =  \.l  as  £  decreases  from  £  =  1  to 
£  =  1/20  (this  corresponds  to  an  increase  of  the  capillary  number). 
Similarly,  for  the  second  phase,  the  power  decreases  from  fi  =  1.4 
to  p  =  1 .3  with  decreasing  £.  On  the  other  hand,  ft  is  much  less  sen¬ 
sitive  to  variations  in  x,  and  remains  almost  constant  for  distinct  x 
(a  change  in  ft  around  0.1 -0.2  for  both  permeabilities  -  results  are 
not  shown  here).  Furthermore,  the  parameter  (B),  which  accounts 
for  the  irreducible  saturation  of  the  first  phase,  approaches  unity  as 
£  decreases.  From  the  fact  that  kr  i  =  0  for  a  phase  (i)  which  does  not 
percolate,  kri  should  be  correlated  as  a  function  of  (s*  -  sim  i)  in  order 
to  obtain  kr,i(si  -  sim,i  =  0)  =  0.  Hence,  the  following  relationship  is 
used: 

kr4  =  Bk{Si-sim^,  (15) 

where  Bk  and  fik  are  again  parameters  to  be  estimated,  and  where 
they  are  influenced  by  x  and  £.  The  shift  in  the  (k^-Si)  curves  from 
Fig.  7  to  account  for  s2m2-  is  shown  in  Fig.  8  (symbols  -  numeri¬ 
cal  results,  and  lines  -  power  law).  The  values  of  the  exponent 
(pk)  decrease  compared  to  j3  and  range  from  1.3  to  1.8  for  the 
first  phase,  while  remaining  essentially  constant  for  the  second 
phase  pk^\.2.  The  relative  permeability  does  not,  however,  fol¬ 
low  the  power  law  over  the  entire  range  of  saturation  differences 
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Fig.  6.  Possible  values  of  (a)  first  phase  relative  permeability  (/<r,i),  (b)  capillary 
pressure  (pc),  and  (c)  first  phase  fluid  front  (A,i)  for  network  random  realizations 
over  entire  range  of  saturations.  These  are  given  with  open  circles;  averages  are 
shown  with  bold  squares  and  corresponding  standard  deviations  with  error  bars 
(horizontal  for  saturation). 


(Sj  —  Simj),  as  can  be  seen  from  the  inset  log-log  plot  in  Fig.  8(a). 
For  small  saturation  differences,  kr  i  deviates  from  the  power  law 
estimated  for  higher  saturation  differences.  As  for  B,  the  values  of 
Bk  are  close  to  one  in  all  cases,  except  for  the  second  phase  and 
f  « 1 ,  suggesting  again  that  for  unstable  flows  the  originally  present 
phase  (the  first  phase),  cannot  be  completely  displaced,  and  hence 

%,i  >0. 

4.3.  Influence  of  x  and  t; 

In  order  to  determine  how  the  (/<ri—  s-pc)  and  (A.,— s)  curves 
change  with  heterogeneity  (/)  and  the  cluster  size  length  scale  (/s) 
(process  extent,  £  =  /s/L),  the  numerical  results  for  all  investigated 
pairs  (/,  f)  are  averaged  (as  shown  in  Figs.  6(a)-(c)  for  x  =  2  and 
£  =  1).  The  relative  permeability  (/<rl)  is  a  phase  permeability  (JQ) 
(i  =  1,  2)  normalized  with  respect  to  the  single-phase  permeability 
(Ksp).  Both  I<i  and  I<sp  are  influenced  by  /.  The  single-phase  perme- 


Fig.  7.  Relative  permeability  as  a  function  of  phase  saturation  for  different  cluster 
size  (symbols):  (a)  first  phase;  (b)  second  phase.  The  lower  and  upper  power  law 
limits  are  plotted  with  dash-dot  and  solid  lines,  respectively. 


Fig.  8.  Relative  permeability  as  a  function  of  phase  saturation  shifted  for  a  phase 
immobile  saturation,  and  lower  (dash-dot)  and  upper  (solid  line)  power  law  approx¬ 
imation:  (a)  first  phase;  (b)  second  phase.  Distinct  cluster  sizes  (£)  are  represented 
by  different  symbols. 
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Fig.  9.  Influence  of  cluster  size  (£)  on  (a)  relative  permeabilities,  (b)  capillary  pres¬ 
sure,  and  (c)  fluid  fronts.  In  part  (c),  the  saturation  range  for  which  both  phases 
percolate  is  given  with  As. 


ability  decreases  as  the  medium  heterogeneity  increases,  whereas 
for  a  drainage  like  process,  where  the  largest  throats  are  occupied, 
I<i  increases  with  x-  There  is  no  influence  of  £  on  I<sp,  but  JQ  depends 
on  f  and  thus,  so  does  krj.  Conversely,  the  capillary  pressure  should 
be  a  strong  function  of  the  porous  medium  heterogeneity  (Eq.  (4)) 
and  of  the  process  extent.  Both  parameters  change  the  throat  size 
distribution  at  the  interface.  For  larger  x,  the  influence  of  small 
throats  is  larger  and  pc  increases;  for  larger  £,  the  phases  are  less 
segregated  and  the  throat  radii  at  the  interface  are  shifted  toward 
larger  values,  and  hence,  to  smaller  pc.  The  fluid  fronts  and  no¬ 
percolation  regime  are  strongly  influenced  by  £  and  to  a  lesser 
degree  by  x- 

Fig.  9  shows  how  all  five  parameters  change  with  saturation 
and  f  for  constant  x  =  2.  The  originally  present  phase  (first  phase) 
and  the  second  phase  are  depicted  by  solid  and  dash-dot  lines, 


respectively.  Different  symbols  are  used  for  distinct  values  of  f. 
When  kri  =  0  for  phase  (z),  there  is  no  percolation,  the  immobile 
part  of  the  phase  does  not  contribute  to  the  momentum  transport, 
and  the  relative  permeabilities  are  given  as  a  function  of  (s,-  -sirn<i) 
(Fig.  9(a)).  As  can  be  observed  from  Fig.  9(a),  there  are  large  dif¬ 
ferences  between  the  two  relative  permeabilities  kr<\  and  kr2  as 
function  of  the  cluster  size  (ls)  (parameter  f).  Both  relative  per¬ 
meabilities  increase  as  f  decreases,  because  the  flow  stabilization 
causes  the  clusters  to  become  smaller.  Hence,  a  larger  part  of  the 
porous  medium  is  available  to  the  flow,  resulting  in  an  increases  in 
kr  i.  This  increase  is  more  prominent  for  the  second  phase,  whereas 
regardless  of  f ,  kTt  i  increases  only  slightly  (the  mobile  clusters  of  the 
first  phase  are  always  present).  The  two  permeabilities,  kr<\  and  fcr>2 
are  not  symmetric,  with  /<r  2  having  maximum  values  smaller  than 
one.  This  is  due  to  the  fact  that  the  second  phase  mobile  clusters  do 
not  exist,  and  the  second  phase  meanders  throughout  the  network 
around  large  clusters  of  the  first  phase  (mobile  or  immobile).  How¬ 
ever,  as  f  decreases  and  the  phases  segregate,  both  phases  behave 
like  the  single  phase,  and  the  relative  permeability  is  influenced 
by  the  porous  medium  volume  occupied  by  that  phase.  Both  rel¬ 
ative  permeabilities  become  more  symmetric,  and  change  almost 
linearly  with  s*. 

The  changes  of  the  capillary  pressure  (pc)  with  f  are  largely  influ¬ 
enced  by  the  heterogeneity  parameter  (x).  Thus,  for  a  network  that 
is  close  to  homogeneous  (all  throats  have  similar  sizes),  the  capil¬ 
lary  pressure  is  constant  with  f.  For  a  heterogeneous  network,  as 
the  preferably  filled  throats  and  the  interface  shape  change  with 
f,  a  variation  in  pc  is  observed  as  shown  in  Fig.  9(b)  for  four  dif¬ 
ferent  values  of  f.  The  capillary  pressure  increases  as  f  decreases, 
where  the  throats  with  smaller  radii  become  part  of  the  interface 
for  smaller  f.  Furthermore,  and  as  observed  earlier,  the  capillary 
pressure  (pc)  does  not  change  with  saturation  when  f  is  held  con¬ 
stant.  In  the  latter  case,  the  only  change  that  occurs  with  varying 
saturation  (s)  is  the  interface  shape,  but  the  average  radius  of  the 
interfacial  throats  remains  constant.  For  small  f ,  the  effects  of  phase 
channeling  are  observed.  Both  phases  percolate,  and  the  fluid  front 
positions  are  Aj  =  l.  As  f  increases,  the  fluid  flow  becomes  unsta¬ 
ble  and  one  phase  can  be  blocked  by  the  other.  The  deficient  phase 
therefore  does  not  reach  the  outlet,  A;  <  1.  Fig.  9(c)  shows  the  fluid 
front  positions  (A;)  where,  for  distinct  f ,  different  saturation  widths 
(As)  are  obtained  for  which  both  phases  percolate.  In  general,  As 
increases  as  f  decreases. 

It  is  instructive  to  examine  how  the  porous  medium  hetero¬ 
geneities  (x)  influence  (kri—  s-pc)  and  (A,,— s)  curves,  and  these 
results  are  shown  in  Fig.  10(a)-(c)  (for  kr  i,  pc,  and  A*,  respectively). 
Again,  as  in  Fig.  9,  the  first  phase  is  plotted  with  solid  lines,  and 
the  second  phase  with  dash-dot  lines.  Different  symbols  are  used 
for  three  values  of  x,  with  the  length  scale  held  constant  f  =  1.  The 
capillary  pressure  is  a  strong  function  of  x-  On  the  other  hand,  a 
comparison  between  kr  l  and  kr  2  reveals  that  kr  \  is  a  weaker  func¬ 
tion  of  x,  whereas  /<r  2  increases  twofold  as  x  changes  from  x  =  0.3 
to  x  =  14-  For  the  first  phase,  this  fact  might  be  attributed  to  the 
mobile  clusters  of  the  first  phase  (originally  present  phase),  where 
the  clusters  are  sufficiently  large  that  they  represent  the  momen¬ 
tum  transport  hindrance  caused  by  heterogeneities  (x)  in  the  same 
manner  as  the  overall  network.  Thus,  one  does  not  observe  changes 
of  kr  i  with  x;  this  is  corroborated  with  percolation  results  for  the 
fluid  front  (Ai)  (solid  lines  in  Fig.  10(c)).  No  noticeable  changes  of 
Ai  with  si  for  different  x  are  observed,  and  the  first  phase  volume 
follows  the  increase  of  the  volume  of  the  overall  network.  There¬ 
fore  s2  =  1  -  Si  is  also  not  influenced  by  x,  and  the  increase  in  kr>2 
may  be  purely  caused  by  the  larger  radii  of  the  second  phase  flow 
paths  with  increasing  x-  The  influence  of  x  on  the  capillary  pres¬ 
sure  is  more  pronounced  than  on  kr  i  and  A.,-.  As  x  increases,  some 
parts  of  the  interface  will  have  smaller  throat  radii  with  larger  capil¬ 
lary  pressures.  Finally,  pc  does  not  depend  on  saturation  for  a  given 
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Fig.  10.  Variation  of  (a)  relative  permeability,  (b)  capillary  pressure,  and  (c)  fluid 
fronts  with  medium  heterogeneity  (x)  for  large  clusters,  £  =  1.  The  second  phase 
relative  permeability,  /<r2  is  for  an  order  of  magnitude  lower  than  kr  d  due  to  the 
mobile  clusters  of  the  first  phase. 


X,  as  interface  radii  average  to  the  same  values  regardless  of  the 
saturation. 

4.3  A.  Displacement  flow 

In  the  limit  when  v  =  ni/n2  becomes  small  (the  second  phase  is 
dominant)  and  approaches  zero,  the  parallel  flow  becomes  a  dis¬ 
placement  flow.  Here,  the  second  phase  is  mobile,  whereas  the 
originally  present  phase  (first  phase)  is  distributed  in  the  network 
as  immobile  clusters.  For  this  case,  the  second  and  first  phase  sat¬ 
urations  are  referred  to  as  the  maximum  saturation  (smaXj2)  and 
irreducible  saturation  (sir  l ),  respectively.  The  saturation  smax>2  is 
found  by  extrapolating  data  from  Fig.  5(a).  Similarly,  the  relative 
permeability  of  the  displacement  flow  (kdiSt  2)  is  obtained  by  extrap¬ 
olating  from  the  (kr>2  ~  s2 )  dependence  in  the  limit  of  s2  approaching 


Fig.  11.  Displacement  flow  relative  permeability  (kdiSi 2):  (a)  variation  with  x  and  £ 
and  power  law  approximations;  (b)  envelope  (kdiSt2-smax,2)  and  the  saturation  power 
laws.  The  saturation  power  laws  differ  from  one  obtained  for  two  mobile  phases 
flow. 


unity.  Thus,  kdis2  becomes  an  intrinsic  value  of  /<r2.  The  perme¬ 
ability  kdiS,2  is  a  function  of  both  x  and  £,  where  it  is  assumed 
kdis,2  =BstJ~Ps  and  Bs  and  /3S  are  parameters  that  depend  on  x*  The 
results  of  such  a  comparison  are  shown  in  Fig.  11(a)  with  extrapo¬ 
lated  kdis2  given  by  the  symbols  connected  with  dashed  lines,  and 
the  power  law  approximation  is  shown  with  solid  lines.  Both  smax>2 
and  the  permeability  kdis2  increase  with  x*  Here,  Bs  and  j3s  depend 
on  x,  and  ps  decreases  as  x  increases.  The  parameter  ps  is  found  to 
be  in  the  range  of  0.7-0.8. 

The  saturation  (smax  2)  can  be  interpreted  as  an  equilibrium  sat¬ 
uration  of  the  displacing  phase.  This  saturation  depends  on  x  and 
£  (Fig.  5(b)),  and  in  order  to  obtain  (kdis2-smaXt2)  curves,  the  satura¬ 
tion  results  shown  in  Fig.  5(a)  are  correlated  with  the  permeability 
data  in  Fig.  11(a).  The  results  are  presented  in  Fig.  11(b),  where  the 
numerical  points  (smaXt2,  kdis  2 )  are  depicted  with  square  symbols, 
and  power  law  fits  are  shown  with  solid  lines  for  constant  x  with 
power  /3X,  and  dash-dot  lines  for  constant  £  with  power  fa.  In  dis¬ 
placement  flow,  fa  is  found  to  vary  in  the  range  of  1 .1  -2.2,  whereas 
the  influence  of  the  cluster  size  f  on  the  power  fa  is  much  higher, 
and  fa  =  4  for  all  x •  The  values  of  the  exponents  imply  that  f  has 
a  larger  influence  than  x  on  the  relative  permeability;  the  same 
trend  was  observed  in  the  flow  of  two  mobile  phases.  It  should  be 
noted  here  that  the  exponents  and  fa  (see  Eq.  (15)),  estimated 
for  the  second  phase  for  flow  of  two  mobile  phases  and  constant  x, 
correspond  to  the  power  fa.  However,  fa  is  larger  than  or  fa,  sug¬ 
gesting  that  different  mechanisms  influence  the  displacement  and 
flow  of  two  mobile  phases.  The  same  relative  permeability  depen¬ 
dence  cannot  be  used  in  both  types  of  flow.  This  distinction  does 
not  apply  to  the  capillary  pressure,  as  pc  does  not  change  with  sat¬ 
uration  in  the  flow  of  two  mobile  phases.  The  envelope  for  pc  is 
therefore  a  function  of  x  and  f  as  shown  in  Fig.  12,  but  not  of  the 
flow  type. 
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Fig.  12.  Changes  of  the  capillary  pressure  (a)  as  a  function  of  medium  heterogeneity 
(x)  and  cluster  size  (£),  and  (b)  construction  of  (pc-s,r,i )  envelope  for  displacement 
flows,  where  sirji  is  irreducible  saturation  of  the  originally  present  phase.  The  enve¬ 
lope  remains  the  same  in  two  mobile  phases  flow  as  pc  is  constant  with  s i  (see  e.g. 
Figs.  9(b)  and  10(b)). 


5.  Conclusions 

Two-phase  flow,  in  which  both  phases  are  mobile  and  flow 
co-currently,  has  been  investigated  using  an  extended  discrete 
pore  network  methodology.  Functional  changes  in  relative  per¬ 
meabilities,  capillary  pressure,  and  fluid  fronts,  as  a  function  of 
phase  saturation,  have  been  predicted  for  conditions  correspond¬ 
ing  to  a  PEMFC  gas  diffusion  layer.  Two  additional  parameters 
have  also  been  varied:  the  network  heterogeneity  and  the  process 
extent,  which  is  a  measure  of  the  maximum  allowed  cluster  size. 
It  was  found  that  the  relative  permeabilities  follow  a  power  law 
of  saturation,  with  the  permeability  of  the  originally  present  phase 
exhibiting  a  higher  sensitivity  to  saturation  changes  than  the  other 
phase.  When  the  immobile  phase  saturation  is  accounted  for  by 
expressing  the  results  in  terms  of  the  reduced  phase  saturation, 
smaller  exponents  are  obtained.  The  numerical  results  suggest  that 
the  cluster  size  influences  the  exponent  much  more  than  the  het¬ 
erogeneity  of  the  porous  medium.  For  a  defined  porous  medium 
and  maximum  cluster  size,  the  capillary  pressure  does  not  change 
with  the  saturation,  but  it  is  found  to  vary  significantly  with  the 


heterogeneity  of  the  medium  and  the  cluster  size.  For  conditions 
that  favor  small  throats  at  the  interface,  the  capillary  pressure 
increases.  This  is  true  for  highly  heterogeneous  media  and  small 
clusters.  The  cluster  size  influences  the  percolation  of  the  phases, 
and  for  small  clusters  the  phase  channeling  occurs.  For  smaller  clus¬ 
ters,  range  of  phase  saturations  for  which  both  phases  percolate 
becomes  broader. 

Displacement  flow  is  a  limiting  case  of  parallel  flow,  where  the 
originally  present  phase  in  the  network  is  highly  deficient  at  the 
network  inlet.  Since  the  numerical  simulations  cover  a  wide  range 
of  saturations,  parallel  flow  data  was  extrapolated  for  the  origi¬ 
nally  present  phase  deficient  at  the  inlet  (second  phase  dominant). 
Relative  permeability  and  capillary  pressure  curves  were  deduced 
from  the  extrapolation.  The  relative  permeability  is  specific  to  the 
flow  type  and  needs  to  be  determined  for  each  flow  type  (parallel  or 
displacement).  In  displacement  flow,  relative  permeability  also  fol¬ 
lows  a  power  law  relationship,  but  the  exponents  are  much  higher 
than  for  parallel  flow.  The  capillary  curves,  on  the  other  hand,  do 
not  depend  on  the  flow  type. 
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